Linear response separation of a solid into atomic constituents: 
Li, Al, and their evolution under pressure 
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We present the first realization of the generalized pseudoatom concept introduced by Ball, and 
adopt the name enatom to minimize confusion. This enatom, which consists of a unique decompo- 
sition of the total charge density (or potential) of any solid into a sum of overlapping atomiclike 
contributions that move rigidly with the nuclei to first order, is calculated using (numerical) linear 
response methods, and is analyzed for both fee Li and Al at pressures of 0, 35, and 50 GPa. These 
two simple fee metals (Li is fee and a good superconductor in the 20-40 GPa range) show different 
physical behaviors under pressure, which reflects the increasing covalency in Li and the lack of it 
in Al. The nonrigid (deformation) parts of the enatom charge and potential have opposite signs 
in Li and Al; they become larger under pressure only in Li. These results establish a method of 
construction of the enatom, whose potential can be used to obtain a real-space understanding of the 
vibrational properties and electron-phonon interaction in solids. 

PACS numbers: 71.15.-m, 71.20.Dg, 74.25. Jb, 74.62.Fj 
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I. INTRODUCTION 

The desire to describe condensed matter as a collec- 
tion of (generalized) atoms is older than the discipline 
of solid state physics itself, extending back to the an- 
cient Greek atomistsiii^ This passion remains as strong 
as ever, embodied in site-based models, atomic orbital 
calculational approaches, and in a resurgence of inter- 
est in Wannier functions. At the other extreme lies the 
homogeneous electron gas starting point for condensed 
systems, where the atomic nature is disregarded at first 
and brought in only as necessary to address the realities 
of real solids. These two viewpoints comprise the tradi- 
tional viewpoints of solid state chemistry and condensed 
matter physics, respectively. The simple model of over- 
lapping free atoms is useful pedagogically but neglects 
the real physics of the formation of solids. 

Anti-intuitively, the itinerant weak-pseudopotcntial 
viewpoint provided the first example of how the den- 
sity n(r) and the electronic potential v(r) can be sepa- 
rated uniquely into atomic contributions. In the limit of 
weak pseudopotentials, the neutral pseudoatom approach 
of Ziman^ gives a description of a nearly-free-electron 
solid as overlapping atomic contributions, a description 
that extends to lattice dynamics and even to the melt. It 
is however severely restricted by the limitation to weak 
pseudopotentials, which applies only to the alkali metals 
and may not give satisfactory accuracy even there. 

Dagens defined an auxiliary neutral atom arising from 
a change in density (with zero net charge) induced by 
a screened potential in jellium. 4 This density functional 
theory inspired prescription was solved numerically for Li 
and Na, but does not address the overlapping of such en- 
tities. In the work of Streitenberger— a generalized pseu- 
doatom model is introduced that extends the pseudoatom 
concept of Ziman to inhomogeneous electron-ion systems 
for simple metals. The model is based on linear-response 



theory in the density functional framework and is applied 
to a metal surface. 



A. Formal background 

Three decades ago Ball introduced a generalized pseu- 
doatom density decomposition concept, one that applies 
to any solid^ It is this specification that we follow in 
this paper. The development of the broad pseudoatom 
concept and the terminology (pseudoatom; auxiliary neu- 
tral atom; generalized pseudoatom; quasi-atom) has a 
long history, and the term pseudoatom also means 'atom 
described by a pseudopotential' and 'phantom atom to tie 
off dangling bonds', as well as many other applications, 
as a literature search will readily reveal. It will there- 
fore be useful to introduce unambiguous language in the 
following: instead of using the generalized pseudoatom 
terminology of Ball we introduce the term enatom^ It 
should be clarified that Ball's pseudoatom has nothing to 
do with any pseudopotential (a common use of the term) . 

For any reference position of atoms, the (vector) first- 
order change in charge density upon displacing one atom 
at Rj from its equilibrium position R®, i.e., the linear 
response to displacement, can be separated into its irro- 
tational and divergenceless components 



dn(r) 
dRi 



-V P3 {r-R°) 



(1) 



VxBjfr-i^); 

here n(r) is the charge density of the system. An imme- 
diate result is a pair of remarkable sum rules^ (i) The 
lattice sum of the rigid density^- pj(r—R°) gives an exact 
decomposition of the crystal charge density into atomic 
contributions: 
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p j (r-R° j )=n(r). 
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(ii) The lattice sum of the deformation density (or "back- 
flow") V x Bj(r — Ft?) vanishes identically: 

J2VxB j (r-R°) = 0. (3) 

3 

This strong constraint reflects that this nonrigid density 
is a cooperative effect of neighboring atoms, which nev- 
ertheless can be broken down uniquely into individual 
atomic contributions. Clearly atoms that are equivalent 
by symmetry have identical pj and Bj ; an elemental solid 
with one atom per primitive cell only has one of each. 

This cooperative origin (a solid state effect) of the de- 
formation density can be understood by considering that 
an atom, embedded in a jellium background, would have 
no deformation density. In fact, if the density is n(r) 
when the atom is at the origin, then by translational 
symmetry the density is n(r — Rj) when the atom is lo- 
cated at Rj. Therefore Vjn(r) = — Vpj(r) and there 
is no nonrigid contribution. Furthermore, within stan- 
dard treatments (such as the local density approxima- 
tion) pj(r) is spherical, because the atom is embedded 
in an isotropic environment. Thus the deformation den- 
sity and the anisotropy of the rigid density arise solely 
from the inhomogeneity of the system, i.e., by neighbor- 
ing atoms. The other extreme is represented by a cova- 
lent solid, which is held together by strong directional 
bonds. The rigid density will be highly non-spherical, 
and the deformation density will be comparatively large, 
reflecting the fact that when an atom moves its bonding 
is disrupted. 

To first order in displacements from the reference point 
5Rj = Rj — R° (typically the equilibrium lattice), the 
density is given by 6 

n(r;{R 3 }) = J2 \Pi ( r " R ° ~ SR i ) ( 4 ) 

3 

+6Rj-Wx Bj(r-R°j)]. 

The quantity inside the sum is the enatom of atom j and 
moves rigidly with the nucleus to first order. (The Rj 

in the argument of Bj can be replaced by Rj without 
changing the expression to first order.) Analogous de- 
composition and sum rules apply to the potential v(r)& 

v(r; {Rj}) = Wi ( r ~ R ° ~ SR 3 ) (5) 

3 

+SRj-S7 x Wj(r-R°j)}. 

Since p, B and V, W are first order quantities, the 
changes in density and potential^ can be related by 
linear response theory^ The deformation arises solely 
from off-diagonal components of the dielectric matrix 
e(q + G,q + G'), i.e., deviation of e(r,r') from the 
e(r — r') form. 

Equation Q allows for a transparent interpretation of 
the quantities Pj and V x Bj. The total charge density 
n(r; {Rj}) of a system of displaced atoms is constructed 



from the charge densities pj(r — R° — SRj) that move 
rigidly with the atoms upon displacement, plus a second 
part V xBj(r—R°) that describes how the charge density 
deforms due to nuclear displacement. 

It is important to keep in mind that, although the rigid 
enatom density (potential) is a specified decomposition of 
the crystal analog, it does not arise simply from screen- 
ing of the pseudopotential (which is the case in weak 
pseudopotential theory). It is intrinsically a dynami- 
cally determined quantity, involving only linear response. 
Specifically, the enatom potential arises from a screened 
displaced (pseudo)potential 

Vj-u = e _1 VjUp S , (6) 

while the enatom density arises from the linear change in 
wave functions 

(occ \ 
2i?e£>*„V^*n • (7) 
kn I 

which can be obtained from first-order perturbation the- 
ory. (The first factor of 2 is for spin.) 

A related quantity is the atomic deformation potential 
Vj£fc n (change in any band energy due to displacement 
of the atom at Rj), given from perturbation theory by 

Vj-e/cn = (kn\S/ jv\kn) (8) 

in terms of enatom quantities. Khan and Allen showed 
the relation of this deformation potential to electron- 
phonon matrix elements £ Resurgent interest in electron- 
phonon coupled superconductivity has led to the sug- 
gestion by Moussa and Cohen that this quantity may 
provide insight into strong coupling. 10 

Although Ball was the first to introduce the enatom de- 
composition and begin to make use of it, the importance 
of Vjn(r) had been recognized earlier. Sham empha- 
sized its essence in the formulation of lattice dynamics, 
related it to the shift in potential by the density response 
function^ and tied its integral to effective charges. Its 
application to ionic insulators was extended by Martin, 
who showed that the enatom dipole and quadrupole mo- 
ments are the fundamental atomic entities that underlie 
piezoelectricity^ 

As powerful as the enatom concept is (see discussion 
below), very little use has been made of it. Falter and 
collaborators have adopted a related quasi-ion idea for 
sublattices of multiatom compounds, and used linear re- 
sponse theory (or models) to evaluate sublattice charges 
for Sii 13 i 14 Ball and Srivastava calculated some aspects of 
the rigid and deformation parts of the density in Ge and 
GaAs from bond-stretch distortions^ No calculation of 
single enatom quantities yet exist for any material. 

B. Motivation 

The enatom concept will be particularly important 
in studying and understanding phonons and electron- 
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TABLE I: Structural and electronic properties of fee Li and 
Al as a function of pressure. Except for the calculated pres- 
sure (P), which is in GPa, all the quantities are expressed in 
atomic units. Vo is the theoretical equilibrium volume, a is 
the fee lattice constant, no is the mean density of electrons in 
10~ 2 el/ag 3 ; N(0) is the density of states at the Fermi level 
in states/(spin Ry atom), Ztf is the Thomas-Fermi screen- 
ing length, Ei? is the Fermi energy (the occupied bandwidth), 
and Icf is the Fermi momentum. a Here a is the experimental 
lattice constant at 35 GPa; our theoretical pressure is 30 GPa. 

phonon coupling, which requires only information aris- 
ing from an infinitesimal displacement of atoms (thus, 
linear response). Current implementations of linear re- 
sponse theory calculate the first-order change in potential 
due to a given phonon, and uses periodicity and Bloch's 
theorem to reduce the calculation to a unit cell (still 
time-consuming). This linear response problem must be 
solved separately for each phonon momentum q. Using 
the enatom concept and linear superposition, it is neces- 
sary only to calculate the enatom density and potential 
once (for each inequivalent atom in the primitive cell) and 
perform elementary integrals necessitating only linearly 
superimposed, overlapping enatom potentials to calcu- 
late the phonon frequencies. (We are concerned here 
only with metals; Ball has shown that insulators with 
long-range potentials require extra considerations. 7 ) The 
electron-phonon matrix elements are even easier, as they 
can be reduced to calculating the matrix elements of the 
enatom potential of each inequivalent atom only. This 
might not be quite as easy as it sounds, because the 
enatom may in some cases have to be calculated out to 
a distance of several shells of neighbors to obtain conver- 
gence of the integrals. We postpone the phonon problem 
to future work. 

For a first detailed application of this enatom con- 
cept to enhance understanding of bonding and electron- 
phonon coupling, we choose the simple metals Li and 
Al. Lithium has attracted renewed interest due to the 
recent discovery that, in spite of being a simple free- 
electron-like metal that is not superconducting above 100 
/iK at ambient pressure^ it displays high T c « 15-17 
K in the 30-40 GPa range,! 7 -^^ and T c = 20 K has 
been reported^ 7 - around 50 GPa. This discovery made Li 
the best superconducting elemental metal (now equaled 
by yttrium 2 ^ and apparently surpassed by calcium 2 ^). 
The evolution of the electronic structure and electron- 
ion scattering within the rigid muffin-tin approximation 
is well studiedj 23 ' Application of microscopic supercon- 
ductivity theory, with phonon frequencies and electron- 



phonon (EP) matrix elements calculated using linear re- 
sponse methodS) 25 ' 26 i 27 ' 28 has established that this re- 
markable level of T c results from strong increase of EP 
coupling under pressure. The one aspect of the electron- 
phonon behavior in Li that is not yet understood^! is the 
strong branch dependence of electron-phonon matrix ele- 
ments. Application of enatom techniques promises to be 
an ideal way to approach the remaining questions. 

Aluminum is the simplest trivalent metal, with T c 
= 1.2 K at ambient pressure; superconductivity is sup- 
pressed with pressure, with T c < 0.1 K at 6 GPa. 29 Un- 
der pressure the electronic structure of Al remains that 
of a free electron-like metal, and a structural transition 
to a hep phase takes place only at P > 217 GPa. 30 Li, 
on the other hand, becomes more and more covalent and 
undergoes several phase transitions^ In the case of met- 
als we use the term "covalency" in a loose sense to in- 
dicate the appearance of directional bonds. The vibra- 
tional and electronic properties of the two systems also 
display important differences. While the electronic struc- 
ture, Fermi surface and vibrational spectrum of Al follow 
a completely normal trend (i.e., the band dispersion be- 
comes steeper, the Fermi surface is virtually unchanged 
and the phonon spectrum is hardened), in fee Li the 
FS evolves from a typical s-like shpere into a multiply- 
connected (Cu-like) shape, with necks extending through 
the L points, reflecting the increase in p character. In the 
phonon spectrum, structural instabilities appear around 
35 GPa along the T — K line, due to the strong e-ph cou- 
pling of some selected phonon modes, whose wave vector 
q connects the necks on the Fermi surface (cf. Refs. [25l . 
\2w . We therefore expect that also the enatom of the two 
systems will display different behaviors under pressure. 



II. CALCULATIONAL APPROACH 

The enatom could be computed by evaluating the 
linear response of a system to the displacement of a 
single atom, that is, the dielectric response e(r,r') <-» 
e(q + G, q + G'). While this may be a viable approach, 
it is demanding and tedious and we use another method 
that requires only minor additional codes. Our approach 
is to let the computer do the linear response for us, by 
choosing a supercell, displacing one atom (taken to be 
at the origin), and obtaining the linear changes Vj-n(r) 
and Vjti(r) by finite differences. Using cubic symmetry, 
displacement in a single direction is sufficient to obtain 
the full vector changes. The enatom components (rigid 
and deformation) are then obtained using the Helmholtz 
decomposition of a vector field into its irrotational and 
divergenceless parts 4 

We obtained the enatom for fee Li and fee Al at atomic 
volumes corresponding to 0, 35, and 50 GPa pressure. 
We used a cube supercell of lattice constant A=3a (lat- 
tice constants are listed in Table H|, which contained 
4x3 3 =108 atoms. The enatom is represented as a Fourier 
series in the supercell. For the jellium calculations we set 
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FIG. 2: (Color online) Evolution of the rigidity factor 1Z 
(Eq. JT0}) of the enatom for Li and Al as a function of pres- 
sure. The difference is the decrease in rigidity with pressure 
in Li, which is magnified somewhat by Li's larger compress- 
ibility. 



FIG. 1: (Color online) Comparison of the magnitude of the 
rigid (top) and deformation parts (bottom) in the first oder 
change of the density (Eq. {l}) of Li at 35 GPa and a jellium 
model. The plot is taken along a [110] direction through the 
atom. The magnitude of |V x B| in jellium is three orders of 
magnitude smaller than in Li, and therefore invisible on the 
scale of the plot. 



a single Li or Al atom into cubic supercells whose lat- 
tice constants correspond to the P=35 GPa cases. The 
mean electronic density was made equal to the related 
crystalline systems and a homogeneous positive jellium 
background provided charge neutrality. 

For the self consistent density functional calculations 
we employed the PWSCF code^ and Troullier-Martins^ 3 . 
norm conserving LDA pseudopotentials and a plane wave 
cutoff energy of 20 Ry for both Li and Al. For the k-space 
integration in the primitive fee unit cell we used a (18) 3 
Monkhorst-Pack grid^ with a cold-smearing parameter 
of 0.04 Ry— With these parameters, we obtained a con- 
vergence of 0.2 mRy in the total energy and of 0.2 GPa 
for the pressure at 35 GPa for both systems. For the 
large cubic supercell, we used a 2 3 Monkhorst-Pack mesh, 
yielding four points in the irreducible Brillouin zone, and 
the same cold-smearing parameter of 0.04 Ry. With this 
choice, the total energy (pressure) calculated in the su- 
percell equals that of the original fee lattice within 0.1 
mRy and 0.1 GPa, respectively. The pressure was calcu- 
lated from a Birch-Murnagham fit of the LDA E vs V 

9? 44 

curve i i 

In Table [J we summarize the most relevant properties 
of Al and Li as a function of pressure. Since the inde- 
pendent variable in our calculations is the volume of the 
unit cell, we also include a column showing the relative 
volume change. We notice that the lattice constant of Li 
decreases very rapidly with pressure, as signaled also by 
the very small bulk modulus. At P = 30 GPa, the unit 
cell volume of Li is already one half of its P — value. 



For comparison, the volume of Al at 50 GPa is 73% of 
its zero pressure value. 



III. ANALYSIS 

As a test of our numerical approach we checked that 
the sum rules of Eqs. ((2|) and ([3|) were almost perfectly 
fulfilled. For the jellium enatom the deformation part 
should be exactly zero, as discussed in Sec. II Al The de- 
formation is not identically zero for our jellium enatom 
due to supercell effects. These effects are however very 
minor, viz. the maximum of the jellium deformation den- 
sity is 500 times smaller than the maximum of the Li 
crystal deformation density. 



A. Deformation vs rigid part 

In Fig. Q] we show the magnitudes of the vector fields 
Vp and VxB along a [110] direction. The rigid parts for 
fee Li and the Li-in-jellium model have the same magni- 
tude and overall shape, while the deformation parts are 
very different. For Li |V X B\ is approximately one order 
of magnitude smaller than |Vp|. This behavior reflects 
a general trend expected in simple metals: the deforma- 
tion part is much smaller than the rigid, nevertheless it 
is quite revealing, as we demonstrate below. 

To quantify the strength of the fields we are considering 
in a more precise way, we define the magnitude M. [A] of 
a scalar or vector field A(r) as the root mean square 



M[A] = ^J^t[A{t)]\ (9) 

where is the volume of the supercell. The relative 
importance of Vp and V x B in Eq. (JXJ) can be quantified 
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FIG. 3: (Color online) A kubic harmonic decomposition of the rigid parts p and V of the enatoms of (a) Li and (b) Al at 35 
GPa. The radial expansion coefficients ft, are defined in Eq. (|llfl . 



by defining a rigidity factor as 



K 



M[V P ] 
M[V x B] ' 



(10) 



This ratio is one measure of how rigidly the enatom den- 
sity (or potential, defined analogously) follows the nu- 
cleus. For the perfect jellium enatom M[V x B] = 
and the rigidity factor would diverge. But because of the 
already mentioned supercell effects, we obtain 1Z ~ 3000 
(1400) for the density and 1Z ~ 2500 (1500) for the po- 
tential of jellium Li (Al), which demonstrates again that 
the supercell effects are indeed small. 

In the actual compounds, for the charge, 1Z has similar 
values at zero pressure but decreases by a factor of 3 be- 
tween and 50 GPa in Li, whereas it is unchanging in Al, 
as shown in Fig. [5] For the potential, 1Z is almost a factor 
of 2 smaller in Li than in Al at zero pressure, and again 
decreases with pressure while that for Al remains con- 
stant. This very different behavior in two simple metals 
is further corroboration that Li is increasing in covalency 
with pressure, while Al is not. The large values of 1Z for 
the potential reflects the fact that the change in potential 
(e.g. due to phonons) is dominated by a rigid part, which 
provides justification for a rigid screened ion or the rigid 
muffin tin potential approximation^ 



B. Rigid part 

1. Lattice harmonic decomposition 

To define the degree of sphericity of the rigid density p 
and potential V these scalar functions can be expanded in 



lattice harmonics of full cubic symmetry, identified with 
angular variation L=0, 4, 6, 8...: 



L 



(11) 



where Kl is the kubic harmonic-- built from spherical 
harmonics of angular momentum L (see Ref. 1381 ). Jl 
are the radial expansion coefficients, and r = \r\. In 
our case the functions Kl are normalized according to 
f(KL) 2 dil = 4-7T, which ensures that the first radial ex- 
pansion term is equivalent to the spherical average, i.e., 
/o(r) = l/(47r)/p(r)dn. 

In Fig. [3] we display the L > radial expansion co- 
efficients at 35 GPa, relative to the (obviously much 
larger) spherical part. The ideal jellium enatom is per- 
fectly spherical and thus contains /o only (see Sec. II Al) . 
The decomposition of our supercell jellium model reveals 
small L > terms, due to small supercell effects. The 
magnitudes of /l for Li and Al are comparable to the 
supercell effects in the jellium enatoms, as the supercell 
boundary is approached. Therefore we conclude that the 
lattice harmonic coefficients f^, fe, /g are meaningful 
only within a radius of ~ A/3- 

While the general characteristics and relative signs of 
the L > terms are quite different in Li and Al, in most 
cases their maximum is only ~ 1% or less of the maxi- 
mum of /o, for both density and potential. The only ex- 
ception is the L — 8 lattice harmonic in the density of Li, 
which is surprisingly large around the nearest neighbor 
distance 4.5 as- This anisotropy reflects the 'coopera- 
tive' influence of neighboring atoms in determining the 
enatom character. The relative size of the L = 8 peak 
grows from 8% to 16% of the maximum of /o from to 50 
GPa. With increasing covalency in Li under pressure not 
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FIG. 4: (Color online) Radial plots of the rigid parts of density and (local, / = 2) potential for (a) Li and (b) Al, both at 35 
GPa. We show the isolated atom, the enatom, the enatom of a jellium model, and, only in the lower panels, the Thomas- Fermi 
potential. The plots show the spherical parts /o of these quantities. On the radial direction the eight nearest neighbors are 
located at (a) 4.53, (b) 4.87 as, the six next nearest neighbors are at (a) 6.41, (b) 6.89 as and the edge of the supercell is at 
(a) 13.59, (b) 14.62 as, illustrating the very small overlap from neighboring supercells. In the inset we show a blow-up of the 
tail region, where the Friedel oscillations of the charge are clearly visible (see also inset in Fig. [5j). 



only / 8 but all non-spherical contributions to the rigid 
density increase; this effect cannot be seen in Al. 

For the enatom potentials, the non-spherical terms are 
small enough in fee Li and Al at all volumes studied that 
the enatom potentials can be considered effectively spher- 
ical, as is the common assumption in simple metals 



2. Screening effects 

Figures [4] and [5] show the rigid parts of density and po- 
tential and their pressure evolution. From linear screen- 
ing theory we know the spherical part of the induced 
change in charge density Aro(r) for a simple metal will 
have long-range (but rapidly decaying) Friedel oscilla- 
tions. Our approach reproduces the long-range oscilla- 
tions according to 

An(r) ~ cos(2fc F r)/r 3 . (12) 

with reasonable agreement (see insets in Figs. |4] and [5|; 
they should be important primarily for describing long- 
range force constants. The corresponding oscillations of 
the enatom potential are not expected to be as impor- 
tant, and we confirm this expectation, as the oscillations 
visible in the inset in the lower panel of Fig. 0Ja) are 
indeed very small. 

In the lower panels of Fig.[4]we also see that the screen- 
ing in the solid causes the enatom potentials to be sig- 
nificantly more short ranged than the atomic potentials. 
Furthermore, the enatom potential is less attractive by 



0.7 to 0.9 Ry in these atoms. The gradients, which de- 
termine electron-phonon matrix elements, do not seem 
to differ greatly. Note that the pseudopotential we have 
used is non-local, and in the plots only the local (I = 2) 
component is shown. The total potential will include 
the £ = 0, 1 nonlocal parts, which are non-vanishing only 
within the core radius (r corc ~ 2 as) and move rigidly. 

To understand this screening better we compare the 
rigid enatom potential V with the Thomas-Fermi poten- 
tial 



V TP (r) = -^—-e- r ^, (13) 



where Q = eZ is the total charge of the valence 
electrons and Ztf is the Thomas-Fermi screening length 
calculated from the mean valence density^ (given in Ta- 
ble Hj. For both systems and all pressures we find good 
agreement for the long range behavior, confirming that 
the system is still dominated by homogeneous electron- 
gas screening. The electronic density of Al is higher than 
the one of Li and therefore the screening is stronger. As 
a consequence the effective potential in Al is more lo- 
calized than the one of Li. The agreement with linear 
screening is better for Al than for Li, further support- 
ing the deviation of Li from the homogeneous electron 
density picture. 
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FIG. 5: (Color online) The pressure evolution of the rigid part of the enatom in (a) Li and (b) Al, plotted along the [110] 
direction. The dotted line represents the difference between the quantities calculated at P=35 GPa and those at P=0. Under 
pressure, the enatom density increases at its peak (at 1.5-2 Ob), and the local (Z = 2) potential becomes less attractive due to 
the increased screening. 



3. The rigid enatom density and potential 

The spherical average fo(r) of the rigid enatom density 
p contains a charge equal to the valence, which is com- 
pared to a valence density for the isolated atom in Fig. SI 
and also with the corresponding enatom in jellium of the 
appropriate density. Note first that, while in a pseudopo- 
tential calculation the density (potential) inside the core 
radius does not have much physical meaning, changes 
within the core radius will still be useful probes of the 
enatom character. 

The enatom density and potential are both more lo- 
calized than in the isolated atom. This difference can be 
ascribed to two effects. 

(i) The density in the tail region is screened in the solid, 
making the effective potential more short ranged and 
causing charge to move inward. As a result the peak 
value around 2 as increases. 

(ii) There is also charge that moves outwards from the 
core region, causing the peak value to increase further 
but also to move outwards. 

The tail is very similar in jellium and in the solid (a 
consequence of similar Thomas-Fermi screening), but in 
the core region and around the maximum the densities 
are different. The jellium enatom density becomes nega- 
tive near the nucleus, with the amount of negative density 
in the core region being about 1% of the valence, for both 
the Li and Al jellium models. Thus this region does not 
contain a significant amount of negative density, but it 
still causes some charge to move away from the core. As 
a result the peak value of the jellium enatom is slightly 



higher and further out than the one of the crystal enatom. 
(Note: in actual supercell calculations there is never any 
negative density, as the negative dip is compensated by 
tails of neighboring atoms.) 



4- Pressure evolution 

The pressure evolution of the rigid enatom density and 
potential are compared on an absolute length scale in 
Fig. Here we can identify aspects of the same two 
effects as described in the preceding section. Under pres- 
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FIG. 6: (Color online) Integrated spherical density Q(r) = 
4tt J r r' 2 f (r')dr' for Li at 0, 35, and 50 GPa. The density is 
shown in Fig. [5] For r > 2 as the shift in density inward with 
pressure is evident, although there is little difference between 
35 and 50 GPa. For r < 2 as a small amount of charge is 
also shifted outward. 
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sure the enhanced electronic screening makes the effective 
potentials more short ranged, and results in the screened 
charge moving inward (effect (i)). For Al the increasing 
pressure causes first a decrease of the extent of rigid den- 
sity around 4 a B (see inset) and second, an increase of the 
peak value around 2 as of about 12% from to 50 GPa. 
For Al the potential change with pressure is negligible. 

Effect (i) (screening) has a stronger influence in Li than 
in Al because, first, the density is lower so screening is 
less, and second, the relative volume change is larger. 
The peak value increases by more than 1/3 from to 50 
GPa. 

But here we also find effect (ii) which causes a small 
amount of charge to move away from the core and leads 
to an outward shift of the peak position in the rigid den- 
sity from 1.9 a B (P=0) to 2.1 a B (50 GPa). The shift 
of charge is indicated more clearly in Fig. [6l where the 
integrated charge is shown for Li. 

The pressure evolution of the enatom potential is char- 
acterized by a decrease in attraction at small r, from -1.78 
Ry (P=0) to -1.60 Ry (50 GPa) in Li. This decrease is 
fairly uniform over the region out to 3.2-3.5 a B beyond 
which it becomes negligible. In Al the decrease in attrac- 
tion between P=0 and P=50 GPa is only 2%. For 35 and 
50 GPa the rigid density of the Li enatom is negative in- 
side 1 ag, but the amount of negative density contained 
within that region is less than 1% compared to the total 
valence charge, and does not even show up in the plot of 
the integrated charge in Fig. O 



C. Deformation part 

In general the vector fields B and W (Eqs. ((4]) and 
©) that describe the deformation parts of density and 
potential have similar morphologies that reflect the sym- 
metry of the lattice. B or W form symmetry related 
"donut" swirls centered at different distances along lines 
connecting the central atom and first (Inn) and sec- 
ond (2nn) nearest neighbors, i.e., the crystal axes (see 
Fig. [7](a) , where only swirls around the 2nn are visible) . 
The swirls associated with the Inn and 2nn have oppo- 
site rotational directions. The derived fields VxBor 
V x W are large at the centers of the swirls of B or W 
(see Fig. [3(c)), and they are primarily directed radially. 

It is informative also to view these fields in planes as 
done in Fig.[7^b) or Fig.[7Ud), where the precise position 
and spatial extent of their features can be judged. It 
can be seen, for example, that the donuts pictured in 
Fig. \Tl,&) are nearly centered on the 2nn Li sites and that 
W is oriented perpendicular to the plane, pointing either 
towards the viewer (vectors visible) or in the opposite 
direction (vectors not visible). As only the fields V x B 
or V x W are involved in the calculation of physical 
properties (see Eqs. ([1]), (Q, and ||5j)), we will focus our 
attention on them. A comparison of the deformation 
parts of density and potential for two different pressures, 
for both Li and Al, is given in Figs. [5] and [5] 



1. Density deformation 

In Li at P=0 the maxima of the charge deformation 
V x B, shown in Fig. [5J are strongly localized around 
the nearest neighbors. Under pressure these maxima 
are pulled inward. The direction of the field determines 
the sign of the charge deformation. For Li in Fig. [8J 
SR-X7 x B (see Eq. g}) is positive for SR || [100], so 
there is a 'charge transfer' from behind the displaced 
atom, to in front of it. Such charge distribution re- 
flects a displacement-induced dipolar moment described 
by the deformation (and which will be screened locally 
in a metal). At ±45° (at the Inn sites, in fact) there is 
a depiction of charge, with a corresponding increase at 
±135° (on the Inn behind the displaced atom). In Al 
the pattern is similar, but the sign is reversed and the 
maxima are nearer the nucleus. These differences will af- 
fect their dynamical properties differently; this influence 
may be significant in Li and is probably not in Al as it 
remains more free-electron-like. 

Under pressure, the magnitude of M. [V x B] in Li in- 
creases quite significantly, being 2.8 x 10~ 5 at P=0 and 
increasing by an order of magnitude at 50 GPa. This 
change, consistent with increased covalency, is the cause 
of the large drop of the rigidity factor in Fig. [21 The pres- 
sure evolution in Al is marginal: M[V x B] is 1.9 x 10 -4 
at P=0 and increase by only 30% at 50 GPa. 



2. Potential deformation 

The potential deformation in Li undergoes a surpris- 
ingly large pressure evolution, reflected in the shape, 
magnitude M., and extent of V x W. M[V x W] is 
4.4 x 10~ 4 at P=0 and increases by over a factor of 4 
by 50 GPa. The contour plot of Fig. [5] shows the change 
from P = to 35 GPa. Starting with a small deformation 
located on the Inn, maxima in the deformation grow in 
substantial regions including the 2nn. 

For Al, V x W has its maxima along the cubic axes, 
and much closer to the nucleus. As for the charge defor- 
mation, V x W has the opposite sign compared to Li, 
and its change with pressure is minor. 

Given the simple shape of the deformation term V x 
W, it is easy to understand its effect on the total change 
in the potential (Eq. ([5])). SR-W x W gives an additional 
dipolar-type contribution, adding to the main change of 
potential Vj-u which has a dipolar form arising from dis- 
placement of the (nearly spherical) rigid potential. 

The pressure evolution of the rigidity factor for the po- 
tential in Li (see Fig. [5]) shows that at 50 GPa the defor- 
mation part contributes about 2% to the total change in 
potential V ' jV (for Al this contribution is negligible). For 
materials with lower rigidity the deformation part might 
give substantial contributions to WjV, large enough to af- 
fect its scattering properties or the strength of electron- 
phonon coupling. 
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FIG. 7: (Color online) The deformation part of the local potential of Li at P = 35 GPa. (Left-hand panels) Three dimension 
isocontour graphs of the magnitude, with arrows indicating the direction; (right-hand panels) contour plots in the (001) plane, 
(a) and (b) show the vector field W , (c) and (d) V x W . The dark green (dark gray) balls represent the position of the central 
atom and the nearest and next nearest neighbors within the supercell, which is displayed as black boundary box. The orange 
(light gray) isocontours in (a) and (c) indicate \W\ = 3.5 x 10~ 3 and |V x W\ = 4.8 x 10~ 3 , respectively. The black arrows 
are field vectors that are located on the isocontours. (b) and (d) indicate the magnitude of the vector fields within a;i/-planes 
that are indicated as black-lined squares in the 3D graphs. Superimposed is a mesh of yellow (light gray) field vectors which 
are located within the plane. 



Additional to the figures shown in this paper we pro- 
vide several color graphs, showing examples of enatom 
quantities for Li in 3D and 2D views, as supplementary 
material for download^ 



IV. DISCUSSION AND SUMMARY 

In this paper we have provided a numerical linear 
response approach, and the first explicit examples, of 



the enatom (the generalized pseudoatom introduced by 
Ball£) density and potential for Li and Al, at pressures of 
0, 35, and 50 GPa. This enatom consists of a rigid and 
a deformation density (and potential). The rigid part 
defines a unique decomposition of the equilibrium den- 
sity (potential) into atomic-like but overlapping contri- 
butions that move rigidly, to first order, with the nuclear 
position. The deformation density (potential) describes 
(again to first order in the displacement) how this charge 
(potential) deforms, and can be viewed as a backflow, 
or (depending on its shape) as a mechanism that trans- 
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(a) Li, P = GPa 



x1 r/ (b) Li, P = 35 GPa 
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FIG. 8: (Color online) 2D graphs (see Fig. [TJl of the deformation part of the charge density VxB for Li and Al for and 35 
GPa. Li undergoes a significant pressure evolution arising in the shape and the magnitude of V x B. Al in turn changes only 
slightly. For the meaning of the symbols see Fig. [7] 



fers charge from one side of the displaced atom to the 
other. The enatom quantities were obtained from super- 
cell finite-difference calculations, demonstrating that this 
approach provides a feasible numerical treatment. 

A rigidity factor 1Z was introduced to quantify the rel- 
ative importance of the rigid and deformation parts of 
the enatom, i.e., characterize how rigidly the charge (or 
potential) moves upon displacement. The rigidity factor 
is expected to be smaller for covalent materials whose 
bonding is strongly direction-dependent, and larger for 
metals that lack such strong bonding^ It has been em- 
phasized recently that Li becomes more covalen t 26 ' 27 un- 
der pressure; the various components of the Li enatom 
have substantiated this trend and provided specific ways 
in which this covalency arises. Aluminum, on the other 



hand, remains quite free-electron like up to 50 GPa. Both 
behaviors are clearly reflected in the pressure evolution 
of the rigidity factor of both density and potential: it 
decreases by a factor of 3-4 in Li but stays almost con- 
stant in Al. Moreover, the rigidity of the potential is 
approximately one order of magnitude bigger than for 
the density. This rigidity supports the picture of a rigid 
potential shift with displacement in both Li and Al. In 
general, the rigidity factor 1Z may become a useful tool 
for quantifying a "generalized covalency" of a system, 
even in the case of metals. 

By kubic harmonic decomposition of the rigid enatom, 
we have shown that in Li and Al the potentials are effec- 
tively spherical, providing support for spherical approx- 
imations in rigid-atom models of electron-phonon cou- 



11 




FIG. 9: (Color online) 2D graphs (see Fig. [TJl of the deformation part of the (local) potential V x W for Li and Al for and 
35 GPa. Plots (a),(b) and (c),(d) share the same color bar, respectively. Lithium undergoes a significant change with pressure, 
while aluminum remains almost unchanged. For the meaning of the symbols see Fig. [7] 



pling. Non-spherical contributions in the rigid density 
become larger as Li becomes more covalent under pres- 
sure, and non-rigid contributions to the potential increase 
somewhat in Li. Changes in aluminum are much smaller. 

The basic features of the spherical part of the rigid 
density and potential can be understood by means of 
linear screening theory (and could be calculated in that 
way). First, the localization of the rigid enatom potential 
is a result of free electron-like (Thomas- Fermi) screening, 
showing that the mean radius of Al is smaller than the 
one of Li. Second, the tails of the rigid densities exhibit 
rapidly decaying Friedel oscillations. 



Another finding is that the rigid enatom density is 
more localized than the density of an isolated atom. This 
is a result of two effects, (i) The density in the tail region 
is screened in the solid, making the effective potential 
more short ranged and causing charge to move inward, 
(ii) There is also charge that moves outwards from the 
core region induced by the potential being less deep than 
in the free atom. The second effect is small compared to 
the first. The same two effects also cause an increase of 
localization of the rigid enatom density when the pressure 
is increased. 

These calculations reveal the manner in which the de- 
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formation part of the enatom reflects the symmetry of 
the lattice and undergoes a pressure evolution, which is 
quite significant in Li but small in Al. The basic morpho- 
logical features of the deformation parts are the same in 
Al and in Li and also in the density and in the potential 
but their position, sign and relative magnitude is differ- 
ent. This behavior confirms the expectation that lattice 
symmetry is paramount in determining the character of 
the deformation of the enatom, at least in nearly free 
electron metals. 

The enatom potential is exactly the quantity that de- 
termines electron-phonon matrix elements, and it is for 
the electron-phonon problem that we foresee the impor- 
tant applications of enatom quantities (see Sec. LA and 
LB). The numerical procedure presented here thereby 
provides a viable means to an improved understanding 
of electron-phonon coupling, based on a real space pic- 



ture. Furthermore, the rigidity factor 1Z could be a useful 
tool to quantify a "generalized covalency" of a system, 
even in the case of metals. Applications to strongly cou- 
pled elemental metals and compounds will be presented 
elsewhere. 



V. ACKNOWLEDGMENTS 

The authors thank O. K. Andersen for ideas and 
encouragement. J.K. gratefully acknowledges support 
from the International Max Planck Research School 
for Advanced Materials (IMPRS-AM). W.E.P. was sup- 
ported partially by National Science Foundation grant 
No. DMR-0421810, and is grateful for support from the 
Alexander von Humboldt Foundation. 



1 Often ascribed to Democritus, ca. 430 B.C. 

2 See, for example, C. Bailey, The Greek Atomists and Epi- 
curus (Clarendon, Oxford, 1928). 

3 J. M. Ziman, Adv. Phys. 13, 89 (1964). 

4 L. Dagens, J. Phys. C 5, 2333 (1972). 

5 P. Streitenberger, phys. stat. sol (b) 116, 179 (1983). 

6 M. A. Ball, J. Phys. C 8, 3328 (1975). 

7 M. A. Ball, J. Phys. C 10, 4921 (1977). 

8 W. E. Pickett, J. Phys. C 12, 1491 (1979). 

9 F. S. Khan and P. B. Allen, Phys. Rev. B 29, 3341 (1984). 

10 J. E. Moussa and M. L. Cohen, Phys. Rev. B 74, 094520 
(2006). 

11 L. J. Sham, Phys. Rev. 188, 1431 (1969). 

12 R. M. Martin, Phys. Rev. B 5, 1607 (1972). 

13 C. Falter, W. Ludwig, A. A. Maradudin, M. Selmke, and 
W. Zierau, Phys. Rev. B 32, 6510 (1985); C Falter, M. 
Selmke, W. Ludwig, and K. Kunc, Phys. Rev. B 32, 6518 
(1985). 

14 C. Falter, H. Rakel, M. Klenner, and W. Ludwig, Phys. 
Rev. B 40, 7727 (1989). 

15 M. A. Ball and G. P. Srivastava, J. Phys.: Condens. Matt. 
4, 1947 (1992); ibid. 5, 2511 (1993). 

16 K. M. Lang, A. Mizel, J. Mortara, E. Hudson, J. Hone, M. 
L. Cohen, A. Zettl, and J. C. Davis, J. Low Temp. Phys. 
114, 445 (1999). 

17 K. Shimizu, H. Ishikawa, D. Takao, T. Yagi, and K. 
Amaya, Nature 419, 597 (2002). 

18 V. V. Struzhkin, M. I. Eremets, W. Can, H.-K. Mao, and 
R. J. Hemley, Science 298, 1213 (2002). 

19 S. Deemyad and J. S. Schilling, Phys. Rev. Lett. 91, 167001 
(2003). 

20 J. J. Hamlin, V. G. Tissen, and J. S. Schilling, Phys. Rev. 
B 73, 094522 (2006). 

21 T. Yabuuchi, T. Matsuoka, Y. Nakamoto, and K. Shimizu, 
J. Phys. Soc. Japan 75, 083703 (2006). 

22 K. Syassen and W. B. Holzapfel, J. Appl. Phys. 49, 4427 
(1978). 

23 L. Shi and D. A. Papaconstantopoulos, Phys. Rev. B 73, 
184516 (2006). 

24 N. E. Christensen and D. L. Novikov, Phys. Rev. B 73, 
224508 (2006). 



G. Profeta, C. Franchini, N. N. Lathiotakis, A. Floris, A. 
Sanna, M. A. L. Marques, M. Luders, S. Massidda, E. K. 
U. Gross, and A. Continenza, Phys. Rev. Lett. 96, 047003 
(2006). 

D. Kasinathan, J. Kunes, A. Lazicki, H. Rosner, C. S. Yoo, 
R. T. Scalettar, and W. E. Pickett, Phys. Rev. Lett. 96, 
047004 (2006). 

D. Kasinathan, K. Koepernik, J. Kunes, H. Rosner, and 
W.E. Pickett, |co~nd-mat/0606537| (to be published). 
S. U. Maheswari, H. Nagara, K. Kusakabe, and N. Suzuki, 
J. Phys. Soc. Japan 74, 3227 (2005). 

D. U. Gubser and A. W. Webb, Phys. Rev. Lett. 35, 104 

(1975) . 

Y. Akahama, M. Nishimura, K. Kinoshita, and H. Kawa- 
mura,Y. Ohishi Phys. Rev. Lett. 96, 045505 (2006). 
M. Hanfland, K. Syassen, N. E. Christensen, and D. L. 
Novikov Nature 408, 174 (2000). 

S. Baroni, A. Dal Corso, S. de Gironcoli, P. Giannozzi, 
C. Cavazzoni, G. Ballabio, S. Scandolo, G. Chiarotti, P. 
Focher, A. Pasquarello, et al. URL |http: / /www. pwscf.org 
N. Troullier, J. L. Martins, Phys. Rev. B, 43, 1993, (1991). 

H. J. Monkhorst, J. D. Pack, Phys. Rev. B, 13, 5188, 

(1976) . 

N. Marzari, D. Vanderbilt, A. De Vita, M. C. Payne, Phys. 
Rev. Lett, 82, 3296, (1999). 

G. D. Gaspari and B. L. Gyorffy, Phys. Rev. Lett. 28, 801 
(1972). 

F. C. von der Lage and H. A. Bethe, Phys. Rev. 71, 612 
(1947). 

K. K. Lehmann and C. Callegan, J. Chem. Phys. 117, 1595 
(2002). 

O.K. Andersen, H.L. Skriver, H. Nohl, B. Johansson, Pure 
and Appl. Chemistry 52, 93 (1979). 

See EPAPS Document No. E-PRBMDO-75-090703 for sev- 
eral color figures showing examples of enatom quantities for 
Li in 3D and 2D views. For more information on EPAPS, 



http://www.aip.org/pubservs/epaps.html 

atom'' 



Greek: "en" denotes within or inside; "atom" denotes in- 
divisible part. Therefore "enatom" connotes the indivisible 
part inside a system. 

Our definition of the rigid part p is predominantly positive 
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[like n(r)] and differs in sign from the convention of Ball.— 
The fields B and W are only defined up to a gauge trans- 
formation: B — > B + V0 leads to the same physical defor- 
mation density for any scalar 'potential' <j>- Nevertheless, 
it makes sense to specify B and W uniquely (up to a con- 
stant) by requiring it to be divergenceless. The Helmholtz 
prescription provides this divergenceless field. 8 
To test the applicability of the pseudopotential method to 
high pressures, we calculated the energy vs. volume rela- 
tion for the two systems and fitted it to a Birch-Murnaghan 
equation, to extract the equilibrium volume (Vb) and the 



bulk modulus at zero pressure (-Bo) and its derivative at 
zero pressure (B' ). For Lithium we obtained : Vo = 127.16 
(a B ) 3 , B = 14.9 GPa, B = 3.33. For Al we obtained: V 
= 105.3 (a B ) 3 , B = 82.7 GPa, B = 4.2, in reasonable 
agreement with the experimental values (Vo = 111.2, Bo 
= 72.7 GPa, B = 4.3).— 

Thomas-Fermi screening length: Itf = l/2(7r/3no) , no 
it the mean electronic valence density in atomic units (see 
Table Q). 



